# -*- coding: utf-8 -*-
"""
Created on Thu Apr 19 09:28:07 2012

@author: Daniele
"""

import numpy as np
import math as mat
import random as rd

import matplotlib.pyplot as plt
from pylab import cm
from numpy import exp
import time


#################################################
# Retorna qual valor deve ser selecionado 
# 
#valor  = 0 minimizacao 
# valor  = 1 maximizacao 
# valor = 2 raizes 
##################################33
#MINIMIZACAO
def criterioselecao(a,b):
    if (a<b):#minimizacao
        return(True)
    else: 
        return(False) 

#MAXIMIZACAO 
'''
    if (a>b):#maximo
        return True
    else : 
        return False 
'''
#BUSCA DE zeros 
'''
    if ( disteucli(a)< disteucli(b) ):
        return True
    else : 
        return False 
''' 



#################################################
# Retorna a norma de dois vetores
#
##################################33

def Dist(ent1,posa,ent2,posb):
    tempdb=0.0
    for j in range(0,D):
        tempdb= tempdb+ ((ent1[posa][j]-ent2[posb][j])**2    )
    return (np.sqrt(tempdb))    

#rho precisa de ajuste para garantir uma boa 
# eliminacao da vizinhanca do otimo encontrado 
def Grasp (entrada,pos,epsilon):
    ssoma=0.0;
    beta = 10000.0
    rho = 1.0 
    result = 0.0
    tempr=0.0
    for i in range (0,len(solucoes)):
        tempr = Dist(entrada,pos,solucoes,i)
        #print tempr , entrada[pos] ,solucoes[i] 
        ssoma = ssoma +(np.exp(-tempr)*Csi(tempr,rho)) # se temmp > rho retorna 0 senao retorna 1 
    result = beta*ssoma
    return (result)

#################################################
# Retorna 1 se entrada menor que rho (epsilon) 
# retorna 0 senao 
#
##################################33

def Csi(aa,bb):
    if (aa <= bb):
        return (1.0)    
    else:
        return (0.0)   

def Funcao (entc,p):
    temp=abs(entc[p][0])+abs(entc[p][1])
    temp2=entc[p][0]**2+entc[p][1]**2    
    f=(-1)*temp*mat.exp(-temp2) 
    
    saida = f + (Grasp(entc,p,epsilon))
   # saida = -(entc[p][0]**2 + entc[p][1]**2)
   #
    #if (f!=saida ):    
    #   print f , saida 
    return (saida)            
    
def minimo(x,D):
    
    minim= x[0,D]
    
    for i in range(HMS):
        if(x[i][D] <minim):
            minim=x[i][D]    
    return(minim)

def maximo(x,D):
    
    maximo= x[0,D]
    
    for i in range(HMS):
        if(x[i][D]>maximo):
            maximo=x[i][D]
    return(maximo)


#################################################
# Retorna a distancia euclidiana do vetor de entrada
#
#
##################################33

def disteucli(entb):
    somat=0.0
    for j in range(0,D):
        somat= somat+ (entb[j]*entb[j]) 
    
    return (np.sqrt(somat))    


# acha pos pior e acha pos melhor , deve me retornar 
# no melhor , o mais proximo de zero 
# no pior , o mais longe de zero 

def achapospior(naoordenado):
    
    pospior = 0
    pior = (naoordenado[0])
    for j in range (HMS):
        if(not(criterioselecao(naoordenado[j],pior))):
            pospior=j
            pior = (naoordenado[j])
    return(pospior)
    
def achaposmelhor(naoordenado):
    
    posmelhor = 0
    melhor = (naoordenado[0])    
    for j in range (HMS):        
        if(criterioselecao(naoordenado[j],melhor)):
            posmelhor=j
            melhor = (naoordenado[j])
    return(posmelhor)    

def HM(N1 , D , HMCR ,PAR,  HMS, limvar , epsilon ):
    
    t=0
    x= np.zeros([HMS,D],float) #armazenando os dados no vetor
    HM=np.zeros([HMS,D+1],float) # x , y e funcao (x,y)
    NOVOH=np.zeros([HMS,D+1],float)
    temp = np.zeros(D,float); # temporario
    for i in range(HMS):    
        for k in range (D): 
            x[i,k] = np.random.uniform(limvar[k][0],limvar[k][1]) #armazenando os dados no vetor
    
    #print x  
    
    for k in range(HMS): 
        for i in range (D):
            HM[k,i]= x[k][i]
        temp=Funcao(x,k)
        HM[k,D]= (temp)
        #print Funcao(x,k)     
        #print HM 
    
    #print "Inicial " 
    #print HM 
    minmax= np.zeros([D,D],float)
    
#minmax [0][0] = minimox  
#minmax [0][1] = maximox 
#minmax [1][0] = minimoy  
#minmax [1][1] = maximoy 
        
    for i in range (D) :        
        minmax[i][0]=minimo(x,i) 
        minmax[i][1]=maximo(x,i)
        
    
    #print minmax
    
    sorteio=0.0
    pospior = 0
    posmelhor = 0
    
    while (t<N1):
        t=t+1
        for i in range(HMS):
            sorteio=np.random.uniform(0,1.0 )
            s_int=rd.randint(0,HMS-1)
            if  (sorteio<=HMCR): 
                if(sorteio<=PAR):
                    for k in range (D):
                        NOVOH[i,k]=x[s_int][k]+sorteio*(minmax[k][1]-minmax[k][0])         #np.std(x[:,k]) #((maximox-minimox))
                    temp=Funcao(NOVOH,i)
                    NOVOH[i,D]=(temp)
                
            else: #bounce back 
                for k in range (D):
                    NOVOH[i,k]=minmax[k][0]+sorteio*(minmax[k][1]-minmax[k][0])
                temp=Funcao(NOVOH,i)
                NOVOH[i,D]=(temp)
            #Achar a posicao do pior  f
            
       # print "HM" 
       # print HM
        pospior=achapospior(HM[:,D]) # retorna o pior valor da funcao aplicada no ponto em HM        
       # print "pospior ", pospior
        
       # print "NOVOHM" 
       # print NOVOH
        posmelhor=achaposmelhor(NOVOH[:,D])# retorna o melhor valor da funcao aplicada no ponto em NOVOHM 
       # print "posmelhor ", posmelhor
       
        # se o melhor NOVOH for melhor que 
        # o pior HM troca
        #print HM[pospior,D]
        #print NOVOH[posmelhor,D]
        if (criterioselecao(NOVOH[posmelhor,D],HM[pospior,D])):
            for k in range (D+1):
                HM[pospior,k]= NOVOH[posmelhor,k]
            
    
    posmelhor=achaposmelhor(HM[:,D])
    #print HM 
    #print posmelhor 
    #for k in range (D):
    #   print HM[posmelhor,k]
    return HM[posmelhor]
    
    
N1=500 # numero de geracoes 
HMS=50# numero de solucoes harmonicas 
D=2 #dimensão do problema
HMCR=0.4 # taxa de valores históricos conhecidos 

PAR=0.8 #  taxa de ajuste dos valores conhecidos 

limvar= np.zeros([D,D],float)

limvar[0][0]=-4;  # x inf 
limvar[0][1]=4;   # x sup
    
limvar[1][0]=-4  # y inf 
limvar[1][1]=4;   # y sup

epsilon = 0.00001; 

num_sol=4
solucoes = [] # vetor das solucoes  [0,0] , [0,1] , [1,0] , [1,1] a segunda casa vai até a dimensao do problema   

start = time.clock()
for j in range (num_sol):
    solucao = HM(N1 , D , HMCR ,PAR,  HMS, limvar , epsilon )
    print " solucao " , solucao
    solucoes.append(solucao)
#olha a ultima geracao e busca uma solucao. Verifica que a solucao nao jah está 
# presente no vetor de solucoes      

end = time.clock()

#print solucoes
print "Tempo de algoritmo = ", end - start, "seconds"   